stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -428.006
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 22 -29.9459 8.1469e-05 0.00165913 0.02854 1 25
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -29.95
## b0 5.36
## b1 1.87
## s 2.71
## m0[1] 38.35
## m0[2] 39.07
## m0[3] 22.51
## m0[4] 28.55
## m0[5] 38.94
## m0[6] 13.39
## m0[7] 29.06
## m0[8] 7.33
## m0[9] 10.49
## m0[10] 22.58
## m0[11] 10.50
## m0[12] 11.12
## m0[13] 32.81
## m0[14] 10.93
## m0[15] 22.21
## m0[16] 8.77
## m0[17] 19.08
## m0[18] 10.43
## m0[19] 42.56
## m0[20] 38.11
## y0[1] 38.88
## y0[2] 43.32
## y0[3] 25.03
## y0[4] 33.24
## y0[5] 42.69
## y0[6] 14.89
## y0[7] 30.94
## y0[8] 6.32
## y0[9] 12.85
## y0[10] 20.92
## y0[11] 12.45
## y0[12] 14.17
## y0[13] 32.38
## y0[14] 12.64
## y0[15] 22.58
## y0[16] 3.00
## y0[17] 20.43
## y0[18] 8.41
## y0[19] 38.89
## y0[20] 36.65
## m1[1] 7.33
## m1[2] 11.25
## m1[3] 15.16
## m1[4] 19.08
## m1[5] 22.99
## m1[6] 26.91
## m1[7] 30.82
## m1[8] 34.73
## m1[9] 38.65
## m1[10] 42.56
## y1[1] 8.33
## y1[2] 12.01
## y1[3] 14.24
## y1[4] 15.40
## y1[5] 26.49
## y1[6] 27.86
## y1[7] 32.41
## y1[8] 38.91
## y1[9] 36.88
## y1[10] 49.14
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -30.47 -30.15 1.23 1.03 -32.82 -29.11 1.01 536 808
## b0 5.33 5.34 1.13 1.12 3.47 7.12 1.01 294 386
## b1 1.87 1.87 0.10 0.10 1.71 2.05 1.00 467 975
## s 3.08 3.02 0.55 0.52 2.32 4.09 1.00 1159 984
## m0[1] 38.34 38.33 1.14 1.10 36.43 40.23 1.00 1852 1199
## m0[2] 39.06 39.05 1.17 1.13 37.11 41.03 1.00 1777 1246
## m0[3] 22.48 22.50 0.69 0.67 21.35 23.59 1.00 750 970
## m0[4] 28.53 28.53 0.78 0.72 27.24 29.82 1.00 1977 1081
## m0[5] 38.93 38.92 1.17 1.12 36.99 40.89 1.00 1790 1230
## m0[6] 13.36 13.36 0.83 0.80 12.03 14.70 1.01 338 491
## m0[7] 29.04 29.04 0.79 0.73 27.73 30.35 1.00 2101 1077
## m0[8] 7.30 7.31 1.05 1.06 5.59 8.98 1.01 296 389
## m0[9] 10.46 10.46 0.92 0.94 9.00 11.94 1.01 313 388
## m0[10] 22.56 22.58 0.69 0.67 21.43 23.67 1.00 759 970
## m0[11] 10.47 10.47 0.92 0.93 9.01 11.95 1.01 313 386
## m0[12] 11.09 11.09 0.90 0.91 9.66 12.54 1.01 319 413
## m0[13] 32.79 32.80 0.91 0.87 31.27 34.31 1.00 2433 1076
## m0[14] 10.90 10.91 0.91 0.91 9.46 12.36 1.01 318 419
## m0[15] 22.18 22.19 0.69 0.67 21.03 23.29 1.00 718 970
## m0[16] 8.74 8.75 0.99 1.00 7.14 10.34 1.01 303 405
## m0[17] 19.06 19.06 0.70 0.67 17.93 20.19 1.00 489 777
## m0[18] 10.40 10.40 0.93 0.94 8.93 11.88 1.01 313 388
## m0[19] 42.55 42.54 1.33 1.30 40.32 44.82 1.00 1475 1295
## m0[20] 38.09 38.09 1.13 1.09 36.20 39.97 1.00 1877 1155
## y0[1] 38.30 38.33 3.30 3.24 32.84 43.56 1.00 1968 1895
## y0[2] 39.01 38.91 3.30 3.03 33.62 44.47 1.00 1591 1758
## y0[3] 22.56 22.58 3.26 3.09 17.22 27.94 1.00 2023 2013
## y0[4] 28.57 28.48 3.26 3.10 23.31 33.81 1.00 1958 2013
## y0[5] 38.91 38.81 3.34 3.16 33.40 44.36 1.00 2056 1904
## y0[6] 13.44 13.49 3.27 3.15 8.00 18.68 1.00 1251 1799
## y0[7] 29.13 29.08 3.20 3.05 23.89 34.35 1.00 1990 1901
## y0[8] 7.38 7.40 3.25 3.03 2.05 12.68 1.00 1473 1680
## y0[9] 10.46 10.50 3.20 3.05 5.30 15.64 1.00 1224 1623
## y0[10] 22.55 22.53 3.28 3.11 17.20 28.16 1.00 1687 2055
## y0[11] 10.50 10.51 3.34 3.38 4.72 15.90 1.00 1549 1808
## y0[12] 11.05 10.89 3.22 3.13 6.15 16.33 1.00 1455 1891
## y0[13] 32.71 32.76 3.20 2.98 27.57 38.02 1.00 2125 1959
## y0[14] 10.86 10.85 3.22 3.12 5.66 16.17 1.00 1519 1924
## y0[15] 22.21 22.16 3.27 3.04 16.92 27.79 1.00 1992 1835
## y0[16] 8.79 8.83 3.30 3.02 3.36 14.25 1.00 1263 1770
## y0[17] 19.12 19.12 3.26 3.12 13.65 24.39 1.00 1942 1974
## y0[18] 10.27 10.39 3.16 2.91 5.07 15.35 1.00 1699 1815
## y0[19] 42.53 42.54 3.36 3.06 37.09 48.02 1.00 1817 1838
## y0[20] 38.04 38.11 3.28 3.20 32.74 43.38 1.00 2046 2009
## m1[1] 7.30 7.31 1.05 1.06 5.59 8.98 1.01 296 389
## m1[2] 11.22 11.22 0.90 0.90 9.80 12.65 1.01 320 424
## m1[3] 15.13 15.13 0.78 0.76 13.88 16.40 1.01 365 563
## m1[4] 19.05 19.05 0.70 0.67 17.92 20.19 1.00 488 777
## m1[5] 22.97 22.99 0.69 0.67 21.83 24.10 1.00 808 1002
## m1[6] 26.88 26.90 0.74 0.69 25.65 28.10 1.00 1551 1057
## m1[7] 30.80 30.81 0.85 0.80 29.41 32.22 1.00 2414 996
## m1[8] 34.72 34.73 0.99 0.94 33.07 36.38 1.00 2265 1127
## m1[9] 38.64 38.63 1.15 1.12 36.71 40.56 1.00 1821 1198
## m1[10] 42.55 42.54 1.33 1.30 40.32 44.82 1.00 1475 1295
## y1[1] 7.33 7.31 3.22 3.10 2.02 12.66 1.00 1518 1746
## y1[2] 11.18 11.20 3.14 3.09 6.10 16.24 1.00 1415 1732
## y1[3] 15.14 15.01 3.12 2.96 10.09 20.35 1.00 1687 1655
## y1[4] 19.12 19.12 3.34 3.12 13.70 24.73 1.00 1965 1892
## y1[5] 23.00 23.11 3.24 3.04 17.81 28.23 1.00 1915 1845
## y1[6] 26.88 26.81 3.19 3.03 21.69 32.12 1.00 1950 1701
## y1[7] 30.81 30.81 3.21 3.00 25.52 36.14 1.00 2220 2033
## y1[8] 34.72 34.75 3.22 3.14 29.61 39.86 1.00 2009 1930
## y1[9] 38.54 38.48 3.38 3.19 32.87 43.77 1.00 2120 1991
## y1[10] 42.35 42.47 3.31 3.18 36.77 47.51 1.00 2067 1925
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 finished in 0.5 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 0.6 seconds.
## Chain 4 finished in 0.5 seconds.
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.8 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.6 seconds.
## Total execution time: 0.8 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -26.37 -28.99 12.60 10.26 -42.98 2.70 1.06 66 34
## b0 5.29 5.32 1.36 1.26 2.89 7.48 1.03 193 37
## b1 1.87 1.87 0.12 0.11 1.67 2.09 1.03 180 33
## s 2.20 2.29 0.97 1.06 0.59 3.72 1.04 79 147
## sx 1.01 1.03 0.60 0.71 0.11 1.99 1.07 38 35
## x0[1] 18.15 17.98 0.85 0.70 16.99 19.78 1.03 243 97
## x0[2] 17.40 17.52 0.86 0.88 15.91 18.56 1.02 264 1265
## x0[3] 8.46 8.56 0.85 0.90 7.01 9.61 1.02 211 755
## x0[4] 12.28 12.35 0.64 0.52 11.17 13.25 1.03 1387 2321
## x0[5] 17.45 17.58 0.82 0.75 16.00 18.64 1.01 425 1486
## x0[6] 3.25 3.42 1.09 1.24 1.25 4.60 1.04 89 142
## x0[7] 13.99 13.77 1.26 1.49 12.43 16.22 1.04 86 251
## x0[8] 0.90 0.98 0.72 0.56 -0.40 2.02 1.02 731 1163
## x0[9] 3.18 3.05 0.80 0.66 2.05 4.60 1.01 517 428
## x0[10] 9.44 9.38 0.70 0.57 8.35 10.60 1.01 820 1865
## x0[11] 1.56 1.75 1.18 1.36 -0.54 3.02 1.04 83 116
## x0[12] 3.98 3.84 0.95 1.10 2.73 5.66 1.02 216 526
## x0[13] 14.84 14.76 0.70 0.56 13.72 16.02 1.02 818 1283
## x0[14] 3.43 3.29 0.78 0.68 2.33 4.78 1.01 421 1578
## x0[15] 9.89 9.76 0.96 1.07 8.63 11.53 1.03 141 558
## x0[16] 1.63 1.72 0.71 0.58 0.38 2.70 1.01 596 1370
## x0[17] 7.98 7.87 0.81 0.85 6.85 9.35 1.02 278 758
## x0[18] 2.34 2.48 0.77 0.64 0.96 3.43 1.02 398 591
## x0[19] 20.01 19.95 0.79 0.58 18.78 21.48 1.05 485 177
## x0[20] 17.31 17.39 0.73 0.61 16.07 18.49 1.02 623 526
## m0[1] 39.18 39.34 1.52 1.48 36.57 41.54 1.00 667 2460
## m0[2] 37.80 37.68 1.77 1.91 35.19 40.71 1.03 87 638
## m0[3] 21.13 21.17 1.60 1.78 18.64 23.61 1.02 278 1790
## m0[4] 28.25 28.30 1.22 1.05 26.20 30.21 1.01 2298 2079
## m0[5] 37.91 37.84 1.69 1.82 35.42 40.71 1.02 137 718
## m0[6] 11.40 11.40 1.96 2.38 8.42 14.36 1.02 223 1917
## m0[7] 31.41 31.25 2.25 2.91 28.17 34.82 1.03 131 1101
## m0[8] 7.00 7.01 1.44 1.29 4.49 9.33 1.02 472 41
## m0[9] 11.25 11.33 1.55 1.56 8.58 13.56 1.02 192 44
## m0[10] 22.93 22.91 1.31 1.17 20.90 25.04 1.01 764 1911
## m0[11] 8.26 8.25 2.11 2.55 5.10 11.58 1.02 216 1389
## m0[12] 12.74 12.79 1.92 2.24 9.45 15.56 1.04 87 36
## m0[13] 33.01 33.05 1.28 1.18 30.88 35.08 1.01 2718 2539
## m0[14] 11.72 11.79 1.57 1.50 8.92 14.16 1.02 155 37
## m0[15] 23.77 23.72 1.78 2.14 21.16 26.46 1.03 115 1771
## m0[16] 8.37 8.36 1.39 1.24 6.12 10.66 1.02 502 115
## m0[17] 20.20 20.21 1.57 1.73 17.62 22.63 1.02 150 91
## m0[18] 9.69 9.64 1.42 1.35 7.53 12.01 1.01 485 1817
## m0[19] 42.66 42.68 1.50 1.37 40.16 45.12 1.01 840 729
## m0[20] 37.63 37.57 1.43 1.37 35.33 40.00 1.01 563 2365
## y0[1] 39.18 39.45 2.77 2.32 34.24 43.30 1.00 2366 3448
## y0[2] 37.83 37.40 3.04 2.69 33.54 43.34 1.01 547 1345
## y0[3] 21.13 20.76 2.83 2.54 17.10 26.12 1.01 875 2472
## y0[4] 28.26 28.20 2.72 2.24 23.97 32.84 1.01 4028 3276
## y0[5] 37.91 37.58 2.97 2.60 33.48 43.28 1.01 675 2143
## y0[6] 11.41 10.96 3.10 2.97 7.17 17.10 1.01 505 2636
## y0[7] 31.36 31.83 3.34 3.28 25.37 35.84 1.02 334 1723
## y0[8] 7.00 6.94 2.82 2.33 2.32 11.64 1.01 1439 698
## y0[9] 11.25 11.51 2.86 2.50 6.22 15.52 1.01 1077 1096
## y0[10] 22.90 23.07 2.70 2.25 18.35 27.17 1.01 4267 3342
## y0[11] 8.27 7.77 3.25 3.11 3.75 14.11 1.01 512 2299
## y0[12] 12.82 13.28 3.04 2.71 7.10 17.08 1.01 358 1744
## y0[13] 32.98 32.99 2.71 2.31 28.51 37.38 1.01 4057 3615
## y0[14] 11.76 12.06 2.85 2.43 6.71 16.03 1.01 896 1180
## y0[15] 23.73 24.17 3.01 2.78 18.24 28.04 1.01 502 2187
## y0[16] 8.37 8.27 2.88 2.30 3.65 13.31 1.01 1775 1309
## y0[17] 20.21 20.53 2.85 2.46 15.16 24.39 1.01 750 2128
## y0[18] 9.68 9.43 2.79 2.37 5.39 14.49 1.01 2317 2746
## y0[19] 42.66 42.65 2.83 2.39 38.01 47.45 1.01 3075 2787
## y0[20] 37.55 37.38 2.76 2.31 33.11 42.31 1.00 2921 2568
## m1[1] 7.26 7.29 1.26 1.17 5.03 9.26 1.03 198 37
## m1[2] 11.18 11.19 1.06 0.99 9.32 12.84 1.02 215 36
## m1[3] 15.10 15.11 0.89 0.86 13.53 16.49 1.02 256 35
## m1[4] 19.02 19.02 0.77 0.76 17.70 20.26 1.01 398 81
## m1[5] 22.93 22.92 0.72 0.70 21.81 24.14 1.01 723 1035
## m1[6] 26.85 26.83 0.76 0.74 25.65 28.10 1.01 872 976
## m1[7] 30.77 30.78 0.87 0.85 29.37 32.18 1.01 672 773
## m1[8] 34.69 34.70 1.03 1.02 32.99 36.31 1.02 457 717
## m1[9] 38.61 38.61 1.23 1.18 36.55 40.61 1.02 343 282
## m1[10] 42.53 42.52 1.44 1.38 40.10 44.85 1.02 274 218
## x1[1] 1.06 1.06 1.17 0.77 -0.90 2.99 1.03 3896 2269
## x1[2] 3.13 3.13 1.20 0.76 1.07 5.12 1.03 3995 2078
## x1[3] 5.25 5.25 1.15 0.79 3.25 7.15 1.03 3794 2527
## x1[4] 7.33 7.34 1.16 0.74 5.41 9.28 1.03 4198 2082
## x1[5] 9.42 9.43 1.19 0.81 7.37 11.39 1.03 3848 2244
## x1[6] 11.52 11.54 1.14 0.75 9.62 13.38 1.03 3596 1975
## x1[7] 13.60 13.62 1.18 0.77 11.61 15.55 1.03 4012 1802
## x1[8] 15.73 15.73 1.17 0.76 13.74 17.71 1.03 3716 1926
## x1[9] 17.83 17.83 1.18 0.80 15.83 19.81 1.03 4280 2634
## x1[10] 19.92 19.92 1.17 0.78 17.94 21.93 1.03 3803 2127
## y1[1] 7.22 7.24 2.69 2.33 2.76 11.57 1.00 1187 655
## y1[2] 11.17 11.16 2.60 2.28 6.94 15.34 1.01 1498 1790
## y1[3] 15.11 15.14 2.59 2.19 10.78 19.36 1.01 2309 2993
## y1[4] 19.06 19.02 2.52 2.09 14.90 23.28 1.01 3240 2807
## y1[5] 22.94 22.91 2.54 2.07 18.85 27.28 1.00 3909 2335
## y1[6] 26.83 26.79 2.61 2.11 22.52 31.14 1.01 3388 2375
## y1[7] 30.71 30.67 2.54 2.12 26.43 34.93 1.01 2912 2913
## y1[8] 34.73 34.68 2.62 2.24 30.39 39.07 1.00 2300 3063
## y1[9] 38.59 38.57 2.65 2.36 34.25 43.03 1.00 1916 2091
## y1[10] 42.44 42.46 2.79 2.52 38.01 47.09 1.01 1003 1584
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -16689
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -33.1588 0.000192288 0.000581333 1 1 47
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -33.16
## b0 5.70
## b1 1.93
## s 3.18
## m0[1] 11.48
## m0[2] 19.72
## m0[3] 9.90
## m0[4] 10.86
## m0[5] 22.74
## m0[6] 15.07
## m0[7] 7.14
## m0[8] 19.06
## m0[9] 7.24
## m0[10] 7.93
## m0[11] 6.06
## m0[12] 6.90
## m0[13] 12.03
## m0[14] 18.47
## m0[15] 9.32
## m0[16] 11.59
## m0[17] 16.29
## m0[18] 7.16
## m0[19] 6.45
## m0[20] 20.68
## y0[1] 15.93
## y0[2] 18.77
## y0[3] 6.73
## y0[4] 9.68
## y0[5] 24.91
## y0[6] 16.19
## y0[7] 8.19
## y0[8] 19.66
## y0[9] 3.80
## y0[10] 11.99
## y0[11] 3.38
## y0[12] 6.69
## y0[13] 10.58
## y0[14] 22.80
## y0[15] 7.16
## y0[16] 10.62
## y0[17] 16.49
## y0[18] 11.78
## y0[19] 10.87
## y0[20] 18.60
## m1[1] 6.06
## m1[2] 7.92
## m1[3] 9.77
## m1[4] 11.62
## m1[5] 13.48
## m1[6] 15.33
## m1[7] 17.18
## m1[8] 19.04
## m1[9] 20.89
## m1[10] 22.74
## y1[1] 8.61
## y1[2] 11.55
## y1[3] 12.65
## y1[4] 9.79
## y1[5] 7.93
## y1[6] 15.34
## y1[7] 19.99
## y1[8] 18.29
## y1[9] 18.94
## y1[10] 27.20
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.60 -33.23 1.33 1.03 -36.20 -32.18 1.01 637 553
## b0 5.71 5.61 1.35 1.28 3.67 8.06 1.02 263 282
## b1 1.93 1.94 0.30 0.29 1.43 2.40 1.01 381 601
## s 3.59 3.49 0.65 0.58 2.73 4.77 1.00 1478 1302
## m0[1] 11.49 11.47 0.86 0.82 10.13 12.94 1.01 427 646
## m0[2] 19.72 19.71 1.40 1.30 17.48 21.99 1.00 1611 1306
## m0[3] 9.91 9.89 0.94 0.89 8.49 11.54 1.02 322 408
## m0[4] 10.87 10.84 0.88 0.84 9.49 12.39 1.01 372 522
## m0[5] 22.75 22.74 1.79 1.73 19.87 25.61 1.00 1058 1189
## m0[6] 15.08 15.06 0.93 0.89 13.58 16.61 1.00 1739 1466
## m0[7] 7.15 7.10 1.19 1.12 5.37 9.24 1.02 269 284
## m0[8] 19.06 19.04 1.32 1.22 16.90 21.16 1.00 1724 1306
## m0[9] 7.25 7.20 1.18 1.12 5.47 9.32 1.02 269 288
## m0[10] 7.95 7.89 1.10 1.05 6.28 9.88 1.02 275 296
## m0[11] 6.08 6.00 1.31 1.25 4.12 8.36 1.02 264 283
## m0[12] 6.91 6.86 1.21 1.14 5.09 9.04 1.02 267 289
## m0[13] 12.04 12.04 0.85 0.81 10.66 13.43 1.01 498 726
## m0[14] 18.47 18.46 1.25 1.16 16.44 20.47 1.00 1828 1361
## m0[15] 9.33 9.31 0.98 0.94 7.85 11.04 1.02 302 373
## m0[16] 11.61 11.59 0.85 0.82 10.23 13.04 1.01 440 667
## m0[17] 16.30 16.29 1.02 0.97 14.61 17.99 1.00 2007 1342
## m0[18] 7.18 7.13 1.18 1.12 5.39 9.27 1.02 269 286
## m0[19] 6.47 6.40 1.26 1.20 4.58 8.67 1.02 265 283
## m0[20] 20.68 20.67 1.52 1.42 18.21 23.14 1.00 1432 1290
## y0[1] 11.53 11.63 3.75 3.63 5.33 17.55 1.00 1518 1719
## y0[2] 19.64 19.70 3.87 3.65 13.23 25.86 1.00 1929 1744
## y0[3] 10.00 10.04 3.81 3.60 3.61 16.18 1.00 1629 1952
## y0[4] 10.77 10.77 3.77 3.58 4.44 16.96 1.00 1615 1847
## y0[5] 22.79 22.57 4.05 3.86 16.32 29.54 1.00 1976 2031
## y0[6] 15.02 15.07 3.84 3.66 8.46 21.16 1.00 1931 1951
## y0[7] 7.04 7.10 3.80 3.72 0.90 13.23 1.00 1387 1587
## y0[8] 19.01 19.02 3.94 3.89 12.53 25.48 1.00 2043 1878
## y0[9] 7.15 7.13 3.78 3.55 1.05 13.44 1.00 1012 1543
## y0[10] 7.86 7.98 3.85 3.63 1.44 13.84 1.00 1512 1893
## y0[11] 5.99 5.98 3.87 3.65 -0.30 12.26 1.00 1413 1824
## y0[12] 6.91 6.95 3.96 3.86 0.49 13.39 1.00 1484 1534
## y0[13] 11.96 11.95 3.69 3.68 5.96 17.82 1.00 1928 1779
## y0[14] 18.30 18.28 3.91 3.69 11.96 24.92 1.00 2026 1719
## y0[15] 9.26 9.30 3.83 3.83 2.94 15.56 1.00 1476 1770
## y0[16] 11.49 11.58 3.76 3.53 5.35 17.34 1.00 2038 1899
## y0[17] 16.32 16.39 3.79 3.65 10.16 22.53 1.00 1979 1922
## y0[18] 7.31 7.26 3.81 3.59 1.07 13.59 1.00 1215 1536
## y0[19] 6.50 6.52 3.89 3.77 0.22 12.71 1.00 1353 1758
## y0[20] 20.59 20.44 3.92 3.84 14.22 27.11 1.00 1923 1772
## m1[1] 6.08 6.00 1.31 1.25 4.12 8.36 1.02 264 283
## m1[2] 7.93 7.88 1.10 1.05 6.26 9.87 1.02 275 296
## m1[3] 9.78 9.76 0.94 0.90 8.35 11.43 1.02 317 393
## m1[4] 11.64 11.62 0.85 0.82 10.27 13.07 1.01 444 667
## m1[5] 13.49 13.48 0.85 0.80 12.11 14.88 1.00 882 1284
## m1[6] 15.34 15.32 0.95 0.91 13.81 16.88 1.00 1821 1422
## m1[7] 17.19 17.19 1.11 1.04 15.35 19.01 1.00 2007 1362
## m1[8] 19.04 19.02 1.31 1.22 16.89 21.13 1.00 1728 1306
## m1[9] 20.90 20.89 1.54 1.45 18.38 23.38 1.00 1373 1258
## m1[10] 22.75 22.74 1.79 1.73 19.87 25.61 1.00 1058 1189
## y1[1] 6.21 6.16 3.82 3.72 0.00 12.68 1.00 1157 1211
## y1[2] 7.83 7.76 3.79 3.74 1.71 14.10 1.00 1471 1908
## y1[3] 9.78 9.78 3.76 3.56 3.65 15.73 1.00 1489 1876
## y1[4] 11.69 11.56 3.74 3.65 5.63 17.78 1.00 1873 1818
## y1[5] 13.56 13.65 3.66 3.46 7.55 19.31 1.00 1792 1835
## y1[6] 15.25 15.18 3.74 3.60 9.17 21.49 1.00 1931 2012
## y1[7] 17.12 17.15 3.80 3.71 10.77 23.03 1.00 1814 1879
## y1[8] 19.10 19.15 3.88 3.61 12.78 25.29 1.00 1807 1972
## y1[9] 20.98 20.95 3.90 3.81 14.45 27.37 1.00 2027 1856
## y1[10] 22.72 22.74 4.02 3.88 16.26 29.27 1.00 1900 1892
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -103.137
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 18 -9.24928 5.46446e-06 5.79417e-05 0.05057 1 33
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -9.25
## b0 4.95
## b1 1.99
## s 0.51
## m0[1] 10.91
## m0[2] 19.41
## m0[3] 9.29
## m0[4] 10.27
## m0[5] 22.52
## m0[6] 14.62
## m0[7] 6.44
## m0[8] 18.73
## m0[9] 6.54
## m0[10] 7.26
## m0[11] 5.33
## m0[12] 6.19
## m0[13] 11.48
## m0[14] 18.12
## m0[15] 8.68
## m0[16] 11.03
## m0[17] 15.87
## m0[18] 6.46
## m0[19] 5.73
## m0[20] 20.39
## y0[1] 12.35
## y0[2] 14.07
## y0[3] 9.38
## y0[4] 10.42
## y0[5] 22.60
## y0[6] 15.03
## y0[7] 6.64
## y0[8] 18.72
## y0[9] 14.66
## y0[10] 6.47
## y0[11] 5.70
## y0[12] 6.02
## y0[13] 11.99
## y0[14] 18.18
## y0[15] 8.40
## y0[16] 11.08
## y0[17] 15.59
## y0[18] 5.81
## y0[19] 8.51
## y0[20] 16.12
## m1[1] 5.33
## m1[2] 7.24
## m1[3] 9.15
## m1[4] 11.06
## m1[5] 12.97
## m1[6] 14.88
## m1[7] 16.79
## m1[8] 18.70
## m1[9] 20.61
## m1[10] 22.52
## y1[1] 0.65
## y1[2] 7.29
## y1[3] 9.15
## y1[4] 12.71
## y1[5] 12.55
## y1[6] 15.05
## y1[7] 17.45
## y1[8] 18.93
## y1[9] 20.80
## y1[10] 22.20
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.50 -11.16 1.35 1.08 -14.17 -10.06 1.02 445 358
## b0 4.95 4.96 0.31 0.31 4.44 5.42 1.01 897 1011
## b1 1.98 1.98 0.07 0.06 1.86 2.08 1.01 1195 1397
## s 0.61 0.60 0.19 0.17 0.36 0.96 1.01 735 284
## m0[1] 10.88 10.88 0.21 0.22 10.52 11.21 1.01 1277 1309
## m0[2] 19.32 19.35 0.34 0.30 18.72 19.85 1.01 2120 1554
## m0[3] 9.26 9.26 0.22 0.23 8.89 9.61 1.01 1063 1208
## m0[4] 10.24 10.25 0.21 0.22 9.88 10.58 1.01 1176 1277
## m0[5] 22.42 22.45 0.43 0.39 21.65 23.10 1.01 1916 1669
## m0[6] 14.56 14.57 0.24 0.22 14.15 14.93 1.00 2090 1526
## m0[7] 6.43 6.44 0.27 0.28 5.97 6.86 1.01 916 1045
## m0[8] 18.64 18.67 0.33 0.29 18.08 19.15 1.01 2162 1555
## m0[9] 6.53 6.54 0.27 0.28 6.08 6.95 1.01 918 1069
## m0[10] 7.24 7.25 0.26 0.26 6.82 7.64 1.01 938 1082
## m0[11] 5.33 5.34 0.30 0.30 4.83 5.79 1.01 899 1012
## m0[12] 6.18 6.19 0.28 0.28 5.72 6.61 1.01 910 1023
## m0[13] 11.44 11.45 0.21 0.21 11.08 11.78 1.01 1386 1363
## m0[14] 18.04 18.05 0.31 0.27 17.50 18.51 1.01 2197 1654
## m0[15] 8.66 8.66 0.23 0.24 8.28 9.03 1.01 1009 1133
## m0[16] 10.99 11.00 0.21 0.22 10.63 11.33 1.01 1297 1327
## m0[17] 15.81 15.82 0.26 0.23 15.36 16.20 1.00 2219 1693
## m0[18] 6.45 6.47 0.27 0.28 6.00 6.88 1.01 917 1045
## m0[19] 5.73 5.74 0.29 0.30 5.25 6.18 1.01 905 1013
## m0[20] 20.30 20.33 0.37 0.33 19.65 20.89 1.01 2050 1574
## y0[1] 11.97 10.89 38.79 0.91 7.39 15.14 1.00 1960 2179
## y0[2] 19.25 19.33 21.11 0.96 15.65 22.87 1.00 1938 1743
## y0[3] 0.26 9.23 401.74 0.96 4.66 13.01 1.00 1887 1764
## y0[4] 9.83 10.23 21.39 0.90 6.13 15.06 1.00 2039 1955
## y0[5] 22.49 22.44 19.31 1.08 18.69 26.65 1.00 2171 1993
## y0[6] 13.80 14.59 123.71 0.95 10.40 18.67 1.00 2069 2103
## y0[7] 5.87 6.45 14.80 0.93 2.60 10.52 1.00 1925 1855
## y0[8] 18.63 18.67 8.19 0.94 15.23 22.71 1.00 2011 2003
## y0[9] 6.49 6.54 11.73 0.90 2.91 10.30 1.00 2171 1973
## y0[10] 7.97 7.21 35.48 0.93 3.27 10.37 1.00 1975 1882
## y0[11] 4.86 5.29 33.48 0.96 1.62 9.55 1.00 1836 1891
## y0[12] 5.65 6.18 42.25 1.00 2.52 10.36 1.00 1932 1934
## y0[13] 10.56 11.46 35.57 0.93 7.60 15.17 1.00 1981 1934
## y0[14] 14.49 18.04 130.19 1.07 14.54 21.84 1.00 1974 1894
## y0[15] 10.97 8.67 83.76 0.99 4.47 12.61 1.00 1949 1932
## y0[16] 11.13 10.99 18.11 0.95 7.06 14.44 1.00 1680 1849
## y0[17] 7.01 15.83 407.45 0.94 12.87 19.29 1.00 2018 2002
## y0[18] 5.85 6.42 12.70 0.98 1.77 9.82 1.00 1848 1892
## y0[19] 5.30 5.71 14.30 0.99 1.80 9.49 1.00 2019 1932
## y0[20] 20.39 20.29 11.99 0.96 16.58 23.90 1.00 1958 1741
## m1[1] 5.33 5.34 0.30 0.30 4.83 5.79 1.01 899 1012
## m1[2] 7.23 7.24 0.26 0.27 6.80 7.63 1.01 937 1082
## m1[3] 9.12 9.13 0.23 0.23 8.75 9.48 1.01 1050 1193
## m1[4] 11.02 11.03 0.21 0.22 10.66 11.36 1.01 1303 1327
## m1[5] 12.92 12.93 0.22 0.21 12.55 13.26 1.00 1742 1618
## m1[6] 14.82 14.83 0.24 0.22 14.41 15.20 1.00 2129 1503
## m1[7] 16.72 16.74 0.28 0.25 16.24 17.15 1.01 2243 1607
## m1[8] 18.62 18.64 0.33 0.29 18.06 19.13 1.01 2164 1555
## m1[9] 20.52 20.54 0.38 0.33 19.85 21.12 1.00 2035 1547
## m1[10] 22.42 22.45 0.43 0.39 21.65 23.10 1.01 1916 1669
## y1[1] 5.04 5.29 11.61 0.98 1.48 9.04 1.00 1774 1933
## y1[2] 7.31 7.25 20.32 0.88 4.03 11.29 1.00 1663 1739
## y1[3] 7.92 9.13 60.53 0.90 5.85 12.36 1.00 1831 1911
## y1[4] 8.73 11.01 88.95 0.92 7.35 14.87 1.00 2065 2015
## y1[5] 12.12 12.92 25.57 0.91 9.23 17.11 1.00 2228 2001
## y1[6] 15.05 14.82 16.99 0.97 10.47 18.59 1.00 2053 1891
## y1[7] 16.53 16.71 13.64 0.92 12.81 20.26 1.00 2026 1894
## y1[8] 17.85 18.62 18.88 0.99 14.62 21.85 1.00 2071 1971
## y1[9] 21.33 20.52 27.86 1.08 16.73 25.07 1.00 2033 1903
## y1[10] 22.09 22.45 25.17 1.09 18.56 26.29 1.00 1993 1951
censored data
y i=1-N, y~N(m,s)
acutualy ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11.stan
objective variable is censored
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitiondata=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1957.23
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 30 -12.1858 0.000516082 1.11285e-06 1 1 42
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -12.19
## b0 12.48
## b1 2.52
## s 2.78
## m0[1] 20.66
## m0[2] 21.02
## m0[3] 18.89
## m0[4] 23.78
## m0[5] 22.97
## m0[6] 30.78
## m0[7] 22.05
## m0[8] 23.67
## y0[1] 20.27
## y0[2] 22.04
## y0[3] 20.26
## y0[4] 24.80
## y0[5] 26.36
## y0[6] 27.19
## y0[7] 19.73
## y0[8] 23.66
## m1[1] 12.60
## m1[2] 15.06
## m1[3] 17.52
## m1[4] 19.98
## m1[5] 22.45
## m1[6] 24.91
## m1[7] 27.37
## m1[8] 29.83
## m1[9] 32.30
## m1[10] 34.76
## y1[1] 15.40
## y1[2] 15.01
## y1[3] 14.48
## y1[4] 18.98
## y1[5] 20.79
## y1[6] 26.78
## y1[7] 22.70
## y1[8] 29.94
## y1[9] 29.30
## y1[10] 35.01mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.14 -12.74 1.55 1.35 -16.09 -11.37 1.01 524 580
## b0 12.34 12.45 5.11 4.73 3.67 20.88 1.03 206 230
## b1 2.57 2.56 1.19 1.10 0.59 4.48 1.03 223 269
## s 4.26 3.86 1.79 1.29 2.36 7.52 1.00 969 1052
## m0[1] 20.66 20.65 1.89 1.68 17.57 23.70 1.02 327 499
## m0[2] 21.03 21.00 1.81 1.60 18.08 23.91 1.01 376 615
## m0[3] 18.86 18.89 2.43 2.17 14.76 22.84 1.02 234 299
## m0[4] 23.84 23.80 1.69 1.48 21.19 26.53 1.00 1205 1053
## m0[5] 23.01 22.98 1.62 1.38 20.45 25.53 1.00 1086 1007
## m0[6] 30.95 31.01 4.14 3.74 24.42 37.56 1.02 324 427
## m0[7] 22.07 22.05 1.65 1.39 19.38 24.66 1.01 649 879
## m0[8] 23.72 23.69 1.68 1.46 21.07 26.34 1.00 1204 1077
## y0[1] 20.59 20.60 4.86 4.09 12.44 28.11 1.01 1061 1568
## y0[2] 21.16 21.06 4.63 4.09 13.99 28.65 1.00 1460 1825
## y0[3] 18.78 18.81 5.23 4.39 10.03 26.86 1.00 809 1304
## y0[4] 23.92 23.89 4.83 4.07 16.40 31.49 1.00 1948 1595
## y0[5] 23.25 23.19 4.92 3.99 15.47 31.34 1.00 1922 1894
## y0[6] 31.01 31.01 6.25 5.39 20.90 40.74 1.01 624 1092
## y0[7] 22.11 22.17 4.92 4.22 14.61 29.60 1.00 1743 1693
## y0[8] 23.68 23.66 4.92 4.01 15.58 31.23 1.00 2221 1543
## m1[1] 12.46 12.56 5.06 4.67 3.86 20.92 1.03 206 224
## m1[2] 14.97 15.06 3.97 3.61 8.03 21.45 1.03 207 227
## m1[3] 17.47 17.51 2.94 2.64 12.47 22.33 1.03 215 274
## m1[4] 19.97 19.97 2.07 1.78 16.60 23.33 1.02 273 399
## m1[5] 22.48 22.43 1.62 1.36 19.89 24.99 1.00 874 919
## m1[6] 24.98 24.94 1.92 1.66 21.88 27.96 1.00 1015 916
## m1[7] 27.49 27.49 2.73 2.42 23.09 31.70 1.01 484 669
## m1[8] 29.99 30.03 3.73 3.38 24.16 35.91 1.02 345 441
## m1[9] 32.49 32.53 4.81 4.37 24.91 40.02 1.02 300 402
## m1[10] 35.00 35.05 5.92 5.45 25.52 44.30 1.02 278 377
## y1[1] 12.41 12.39 6.93 6.03 1.73 23.03 1.02 343 684
## y1[2] 14.92 14.92 6.22 5.21 4.98 24.95 1.01 389 866
## y1[3] 17.62 17.52 5.72 4.70 9.03 26.87 1.01 510 792
## y1[4] 20.01 19.89 5.25 4.25 11.62 28.00 1.00 1043 1541
## y1[5] 22.69 22.56 4.96 4.23 15.07 30.68 1.00 1586 1770
## y1[6] 24.72 24.82 4.87 4.01 17.00 32.45 1.00 1877 1563
## y1[7] 27.31 27.25 5.37 4.49 19.13 35.86 1.00 1276 1398
## y1[8] 30.00 29.91 5.95 5.05 20.62 39.30 1.01 660 1403
## y1[9] 32.48 32.42 6.78 5.65 21.31 43.36 1.01 474 846
## y1[10] 35.04 35.19 7.50 6.36 22.83 47.20 1.01 429 857
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -177.751
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 27 -19.1875 0.000483383 0.000105904 0.9521 0.9521 33
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -19.19
## b0 5.85
## b1 4.16
## s 4.27
## m0[1] 19.32
## m0[2] 19.92
## m0[3] 16.40
## m0[4] 24.47
## m0[5] 23.13
## m0[6] 35.99
## m0[7] 21.61
## m0[8] 24.28
## y0[1] 16.55
## y0[2] 28.06
## y0[3] 9.68
## y0[4] 27.35
## y0[5] 26.15
## y0[6] 32.04
## y0[7] 18.23
## y0[8] 22.68
## m1[1] 6.04
## m1[2] 10.09
## m1[3] 14.15
## m1[4] 18.21
## m1[5] 22.27
## m1[6] 26.32
## m1[7] 30.38
## m1[8] 34.44
## m1[9] 38.49
## m1[10] 42.55
## y1[1] 7.09
## y1[2] 8.03
## y1[3] 14.02
## y1[4] 20.92
## y1[5] 22.69
## y1[6] 29.69
## y1[7] 28.38
## y1[8] 33.89
## y1[9] 36.67
## y1[10] 41.44mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -19.63 -19.23 1.55 1.24 -22.78 -17.98 1.01 295 264
## b0 1.58 2.59 7.05 6.08 -11.64 10.52 1.04 182 146
## b1 5.23 4.92 1.63 1.32 3.24 8.29 1.03 197 149
## s 6.55 5.82 2.92 2.06 3.57 11.86 1.01 368 345
## m0[1] 18.54 18.77 2.52 2.18 14.23 22.13 1.03 282 225
## m0[2] 19.28 19.47 2.40 2.09 15.26 22.82 1.03 314 250
## m0[3] 14.86 15.25 3.32 2.80 8.86 19.42 1.04 202 175
## m0[4] 25.01 24.81 2.13 1.72 22.09 28.57 1.01 1076 796
## m0[5] 23.33 23.27 2.06 1.74 20.29 26.83 1.01 749 743
## m0[6] 39.51 38.57 5.47 4.28 32.96 49.49 1.01 274 270
## m0[7] 21.41 21.45 2.13 1.82 18.06 24.87 1.02 459 430
## m0[8] 24.77 24.58 2.11 1.70 21.84 28.33 1.01 1040 751
## y0[1] 18.55 18.54 7.75 6.20 6.89 29.77 1.00 1778 1370
## y0[2] 19.36 19.63 7.43 6.41 7.37 30.57 1.00 1275 1216
## y0[3] 14.77 15.24 7.68 6.58 2.26 25.87 1.01 965 872
## y0[4] 24.92 24.84 7.66 6.03 12.68 36.96 1.00 1892 1523
## y0[5] 23.26 23.29 7.27 6.06 11.81 34.60 1.00 1783 1205
## y0[6] 39.66 38.58 9.22 7.34 27.50 55.61 1.01 624 469
## y0[7] 21.08 21.36 7.70 6.10 9.55 32.49 1.00 1564 1097
## y0[8] 24.89 24.66 7.54 6.15 13.55 36.76 1.00 2049 1618
## m1[1] 1.82 2.83 6.98 6.02 -11.28 10.68 1.04 182 146
## m1[2] 6.93 7.71 5.48 4.60 -3.22 13.99 1.04 184 152
## m1[3] 12.03 12.53 4.05 3.46 4.70 17.37 1.04 189 164
## m1[4] 17.13 17.42 2.80 2.37 12.26 21.05 1.04 232 208
## m1[5] 22.24 22.22 2.08 1.77 19.06 25.71 1.02 556 612
## m1[6] 27.34 27.05 2.42 1.93 24.21 31.55 1.01 910 695
## m1[7] 32.44 31.91 3.53 2.79 28.06 38.91 1.01 397 323
## m1[8] 37.55 36.76 4.91 3.85 31.62 46.59 1.01 292 272
## m1[9] 42.65 41.67 6.39 5.05 34.92 54.45 1.02 256 276
## m1[10] 47.75 46.50 7.92 6.35 38.15 62.62 1.02 239 233
## y1[1] 2.01 3.23 9.52 8.07 -14.36 15.15 1.02 386 349
## y1[2] 6.71 7.86 9.19 7.78 -9.60 18.74 1.01 469 405
## y1[3] 12.02 12.48 8.10 6.74 -2.01 23.75 1.01 624 1014
## y1[4] 17.05 17.48 8.00 6.21 3.27 29.08 1.01 871 775
## y1[5] 22.13 22.20 7.48 6.11 9.89 34.17 1.00 1842 1544
## y1[6] 27.43 26.97 7.51 6.23 16.14 40.22 1.00 1730 1061
## y1[7] 32.48 32.04 7.85 6.32 20.76 45.86 1.00 1370 722
## y1[8] 37.75 36.76 8.87 6.96 25.99 52.87 1.00 936 708
## y1[9] 43.00 41.74 9.45 7.58 30.54 60.08 1.00 527 527
## y1[10] 47.63 46.20 10.60 8.59 34.35 66.00 1.01 419 507
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
sensitivity/specificity analysis
ex14.stan
estimating sens and spec
data {
int N;
array[N] int x;
array[N] int y;
}
parameters {
real<lower=0,upper=1> p;
real<lower=0,upper=1> se;
real<lower=0,upper=1> sp;
}
model {
p~uniform(0,1);
se~uniform(0,1);
sp~uniform(0,1);
for (i in 1:N) {
y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
}
}
generated quantities {
real ppv;
real npv;
ppv=se*p/((se*p)+((1-p)*(1-sp)));
npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}n=20
data=list(N=n,
x=sample(0:1,n,replace=T),
y=sample(0:1,n,replace=T))
mdl=cmdstan_model('./ex14.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -24.7567
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -13.4602 0.00123284 6.67477e-06 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.46
## p 0.32
## se 0.60
## sp 0.40
## ppv 0.32
## npv 0.68system.time({
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,ch=F)
})## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -19.41 -19.11 1.34 1.14 -22.09 -17.89 1.00 672 917
## p 0.51 0.50 0.30 0.39 0.04 0.96 1.00 543 386
## se 0.58 0.58 0.14 0.15 0.34 0.80 1.00 2833 1602
## sp 0.42 0.41 0.14 0.15 0.20 0.65 1.00 2813 1509
## ppv 0.51 0.51 0.30 0.39 0.05 0.96 1.01 567 410
## npv 0.49 0.49 0.30 0.41 0.04 0.95 1.00 579 448
## ユーザ システム 経過
## 0.557 0.116 0.695appendix: bimodal distribution, mixed normal distribution
stan
data {
int<lower=0> N;
real y[N];
}
parameters {
real<lower=0, upper=1> theta; //ratio of mix
real mu1;
real mu2;
real<lower=0> sigma1;
real<lower=0> sigma2;
}
model {
for (n in 1:N) {
target += log_mix(theta,
normal_lpdf(y[n] | mu1, sigma1),
normal_lpdf(y[n] | mu2, sigma2));
}
}EM algorithm
y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 3 components:
##
## log-likelihood n df BIC ICL
## -250 100 6 -528 -530
##
## Clustering table:
## 1 2 3
## 33 25 42rst$parameters## $pro
## [1] 0.329 0.254 0.417
##
## $mean
## 1 2 3
## -5.336 0.201 5.245
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 3
##
## $variance$sigmasq
## [1] 1.05plot(rst)
